%This file generates a case of a 20x20 Cov-matrix with up to 10 sub blocks.
%sampled sequences with different levels of noice are generated and used 
%as estimated subblocks. The performance is measured for different number 
%of blocks and different noise levels. Sample length is fixed to 5000.

clear, close all
NumOfStocks = 20;

load data_10_blocks.mat

TestBlocks=[3 5 8 10];
%TestBlocks=[3:10];

error=zeros(length(TestBlocks),1);

[EstVal,EstLimits]=GetBlockEstimates(CovBlocks, Combinations, NumOfStocks);

for iterBlocks=1:length(TestBlocks)

    NumOfBlocks=TestBlocks(iterBlocks);
    
    %Define weights (dummy)
    w=ones(1,NumOfBlocks);
    
    
    %Generate the optimal estimate
     cvx_begin
        %cvx_quiet(true); 
        variable Sigma_hat(NumOfStocks, NumOfStocks) symmetric;

        %Define objective function
        f = 0;
        s_max = 0;
        for q=1:NumOfBlocks
            f = f + w(1,q)*norm(Sigma_hat(Combinations{q}, Combinations{q})-CovBlocks{q},'fro');           
        end
        minimize (f)

        subject to

        Sigma_hat == semidefinite(NumOfStocks)
        
        %Add constraint which assumes that the estimates are within 20% correct
        %and that the signs are correct of the crosscovariances, which is
        %needed to determine a constraint for each element
        posIndex=find(EstVal > 0);
        negIndex=find(EstVal < 0);
        
        (Sigma_hat(posIndex)) < (1.2)*EstLimits(posIndex);
        (Sigma_hat(posIndex)) > 0
        
        (Sigma_hat(negIndex)) > -(1.2)*EstLimits(negIndex);
        (Sigma_hat(negIndex)) < 0
   
        
    cvx_end
    
    %Compare estimate to original matrix
    error(iterBlocks)=norm(Sigma-Sigma_hat,'fro');
    DiffMatrixRel(:,:,iterBlocks)=abs((Sigma_hat-Sigma)./Sigma);
    DiffMatrixAbs(:,:,iterBlocks)=abs((Sigma_hat-Sigma));
end

figure(1)
for i=1:length(TestBlocks)
subplot(2,2,i)
%clims=[0 0.12*max(max(max(abs(DiffMatrixRel(:,:,:)))))];
clims=[0,1.5]
imagesc(DiffMatrixRel(:,:,i),clims);
colormap(gray)
colorbar
title(['Relative estimation error with ', num2str(TestBlocks(i)), ' sub blocks']) 

end

figure(3)
semilogy(TestBlocks,error,'-.k','linewidth',3)
title('Estimation error')
xlabel('Number of subblocks used')
ylabel('norm(S-Shat)')
